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quantitative test of this model, we first show that the simulated fluctuation spectrum of a one- 
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I. INTRODUCTION 

Dendrites are intricate growth patterns that make up 
the microstructure of many important commercial alloys 
[2,||. They develop a complex shape due to the emis- 
sion of secondary branches behind the growing tips of 
primary branches |g|]. A major advance in understanding 
this dynamical process came historically from the insight 
[Q,g| that results of Zel'dovich et al. g] on the stability of 
flame fronts could be extended to other interfacial pat- 
tern forming systems such as dendrites, viscous fingers, 
etc. For dendrites, further developments along this line 
[R-0| led to a physical picture where small noisy pertur- 
bations, localized initially at the tip, become amplified 
to a macroscale along the sides of steady-state needle 
crystals, thereby giving birth to sidebranches ^-O in 
qualitative agreement with some experiments [H. 

This sidebranching mechanism requires some continu- 
ous source of noise at the tip. Therefore, thermal noise, 
originating from microscopic scale fluctuations inherent 
in bulk matter, is the most natural and quantifiable can- 
didate to consider. Langer ||I9 analyzed the amplifica- 
tion of thermal noise along the sides of an axisymmetric 
paraboloid of revolution and concluded from a rough es- 
timate that it is probably not strong enough to explain 
experimental observations, i.e. sidebranches form closer 
to the tip in experiment than predicted on the basis of 
thermal noise amplification. More recently, Brener and 
Temkin |]ll| made the interesting observation that noise 
is amplified faster along the more gently sloping sides of 
anisotropic (non-axisymmetric) needle crystals, leading 
to the conclusion that thermal noise has about the right 
magnitude to fit experimental data. 

There remain, however, several sources of uncertain- 
ties regarding this conclusion. Firstly, calculations of 
noise amplification have been based on a WKB (Wentzel- 



Kramers-Brillouin) approximation which has only been 
tested by comparison ||] with numerical simulations ||] 
for a fixed frequency perturbation localized at the tip. 
Thermal noise is more difficult to analyze because it in- 
volves a wide range of frequencies and is spatially dis- 
tributed. Consequently, current estimates of the side- 
branching amplitude |lfl,nl| involve some overall prefac- 
tor which is only known approximately. Secondly, the 
predicted sidebranching amplitude depends sensitively 
on the non-axisymmetric tip shape which seems to vary 
from system to system. Bisang and Bilgram |13] have 
found that the tip of Xenon dendrites is well fitted by 
the power law x ^ z'^/^ | p^ , p^ (as opposed to x '^ z^/^ 
for a paraboloid), where z is the distance behind the tip 
and X is the radial distance from the growth axis to the 
interface. In contrast, LaCombe et al. jlGJ, find that 
the tip shape of succinonitrile dendrites is well described 
up to 10 p (where p is the tip radius) behind the tip by 
a weak four- fold deviation from a paraboloid x ~ z^'"^ . 
Since the z^/^ power law should only strictly hold far 
behind the tip ||l4| , [l5| , the proposal [Q that it can be 
used to predict the sidebranching amplitude remains to 
be validated beyond the experiments of Bisang and Bil- 
gram [|3[ . Lastly, analyses of sidebranching have so far 
been constrained to a linear regime. Therefore, there re- 
mains the possibility that nonlinearities produce a noisy 
limit cycle where sidebranches drive tip oscillations. 

At present, it appears to be difficult to make further 
progress on these issues without some reliable computa- 
tional approach to accurately simulate dendritic growth 
with thermal noise. Numerical simulations of dendritic 
growth using a phase-field approach are consistent with 
a noise amplification scenario in that sidebranches are 
absent in purely deterministic simulations where the dif- 
fuse interface region is well resolved |l^-p3]. Moreover, 
in certain simulations sidebranching has been induced 



by randomly driving the tip 1 17 Iq| in a fashion which is 
adequate to produce dendritic microstructures, but not 
to investigate quantitatively the physical origin of side- 
branching. In addition, in front-tracking simulations p3[ , 
sidebranching appears to be due to the amplification of 
numerical noise which is difficult to control. 

The first goal of this paper is to demonstrate that the 
phase-field approach |^,^ can be successfully extended 
to study the effect of thermal noise quantitatively. The 
second goal is to use this approach to carry out a quan- 
titative study of sidebranching in order to test the pre- 
dictions of the linear WKB theory of noise amplification 
|^,Q. Here, simulations are restricted to two dimen- 
sions in order to carry out this comparison in the simplest 
non-trivial test case. There are two main reasons to elect 
a phase-field approach to study thermal noise. Firstly, 
this approach has proven extremely successful to simulate 
dendritic growth |17|-pi[. By reformulating the asymp- 
totic analysis of the phase-field model, it has recently 
been possible to lower the accessible range of undercool- 
ing as well as to choose an arbitrary small interface ki- 
netic coefficient [^. In addition, adaptive mesh refine- 
ment methods, used in combination with the reformu- 
lated asymptotics, have pushed the limit of undercooling 
even further towards the experimental range p2| . Sec- 
ondly, the phase-field approach provides a natural frame- 
work to incorporate thermal noise since it is adapted 
from phenomenological continuum models of second or- 
der phase transitions used to study fiuctuations near a 
critical point pffl. Therefore, the formalism to incorpo- 
rate noise into such models already exists. The extension 
to the phase-field model mainly requires the use of the 
fluctuation-dissipation theorem together with an appro- 
priate scaling of parameters to relate the magnitude of 
the noise in the model with the noise that is present in 
an experiment. This straightforward exercise is carried 
out here. An additional issue is the numerical resolution 
of a small amplitude noise which could be masked by the 
numerical noise and/or discretization artifacts that are 
present in simulations. This problem is absent in stud- 
ies of phase-transitions where the bare magnitude of the 
noise is not important. Here, however, this magnitude 
plays a crucial role. Fortunately, we shall find that it is 
possible to resolve accurately a small amplitude noise, of 
magnitude comparable to experiment, provided that the 
spatially diffuse interface region is well resolved. 

In the context of this study, we are naturally led to 
revisit the issue of the relative importance of the noises 
acting in the bulk and at the interface, which was previ- 
ously considered in the context of a sharp- interface model 
pjj . Microscopically, the bulk noise originates from fluc- 
tuations in the heat current in the solid and liquid phases, 
whereas the interface noise originates from the exchange 
of atoms between the two phases (i.e. the attachment and 
detachment of atoms at the interface) . In ref . p^ , it was 
shown by a direct calculation of the equilibrium fluctua- 



tion spectrum of a flat interface that the bulk and inter- 
face noises drive, respectively, long- wavelength (A > A*) 
and short- wavelength (A < A* ) regions of this spectrum, 
where the cross-over length A* = AircD/fiL. Here, c is the 
specific heat per unit volume, D the thermal diffusivity, 
L the latent heat of melting per unit volume, and /i the 
interface kinetic coefficient. On this basis, it was roughly 
estimated that the bulk noise should predominantly drive 
sidebranching whenever A* < As, where, A5 ~ yUJdoJV, 
is the stability length below which perturbations of the 
interface are stable, V is the tip velocity, and do is the 
capillary length. This condition is actually satisfied for 
growth at low velocity where simple estimations allow 
one to conclude that Ag » A* for materials with reason- 
ably fast attachment kinetics. In the phase-field model, 
the bulk and interface noises are represented by Langevin 
forces added to the evolution equations for the tempera- 
ture and phase fields, respectively. It it therefore possible 
to probe the relative importance of these two noises. In 
this paper, we focus on a low velocity limit where the 
bulk noise should be dominant according to the above 
estimate. We observe indeed hat sidebranching is un- 
affected when the noise is switched off in the evolution 
equation for the phase-field, and only the conserved noise 
is kept in the diffusion equation. 

This paper is organized as follows. In Section II, 
we review the sharp-interface equations of solidification 
with thermal noise and certain useful results of fiuctua- 
tion theory. In section HI, we introduce the phase-field 
model and analyze its equilibrium fiuctuation properties, 
which allows us to relate the parameters of this model to 
the known material parameters that enter in the sharp- 
interface model. In section IV we then discuss the numer- 
ical implementation of the model and present the results 
of a detailed numerical test based on comparing the sim- 
ulated and analytically predicted fiuctuation spectra of 
a stationary interface in thermal equilibrium. Next, in 
section V, we present the results of the simulations of 
dendritic growth and a quantitative comparison of the 
sidebranching characteristics (amplitude and sidebranch 
spacing) of a steady-state growing dendrite to the analyt- 
ical predictions of the WKB theory. Finally, concluding 
remarks are presented in section VI. 



II. SHARP-INTERFACE MODEL 

We consider the standard symmetric model with equal 
thermal diffusivities in the solid and liquid phases. The 
incorporation of fluctuations in this model, with reference 
to earlier works, is discussed in detail in Ref. |g^ and we 
only review here the main results. The basic equations 
of the model are given by: 

dtT^ DV^T-V ■] (1) 

LVn^-cDh-(vT\i-VTU)j + cfi-(jlz-Ju) (2) 



I = Tm -Tk 



V 



(3) 



where T{r, t) is the temperature field defined in terms of 
the three-dimensional position vector r = xx+yy+zz, Tj 
is the interface temperature, T^v/ is the melting tempera- 
ture, r = ^TmI L is the Gibbs-Thomson coefficient where 
7 is the surface energy, Vn is the normal velocity of the 
interface, VT|; (VT|s) is the temperature gradient eval- 
uated on the liquid (solid) side of the interface, k is the 
interface curvature, and other parameters were defined in 
section I. The conserved noise, j = jxX+jyjj+jzZ, repre- 
sents the fluctuating part of the heat current, where the 
components, jVn, with (m — x, y, z), are random variables 
uncorrelated in space and time that obey a Gaussian dis- 
tribution. The variance of this distribution 



<Jm(r,t)j„(f^,t') >=2 



DkBT{r,tf 
c 

X dmnS{r^ 7^)6(1^1'), (4) 



is fixed by the requirement that the diffusion equa- 
tion driven by this noise produces, in equilibrium, the 
known distribution of temperature fluctuations in the 
solid and liquid phases, which is a simple application of 
the fluctuation-dissipation theorem. According to basic 
principles of statistical physics [^ , the mean square fluc- 
tuation of the temperature in a small volume AV^ of solid 
or liquid is given by |E8| : 



< AT^ >= 



cAV 



(5) 



where ks is the Boltzmann constant, which is precisely 
the result that one obtains from a simple calculation of < 
AT^ > using Eq. |l| with J defined by Eq. |. Note that, in 
a non-equilibrium situation, the temperature variation in 
the liquid is small compared to the melting temperature, 
such that T{r,t) can be replaced by Tm on the r.h.s. of 
Eq. |. 

Next, to write down the correlation of the non- 
conserved noise that enters in the interface condition 
(^, it is convenient to define the interface position, 
C(ri,i) = Zint{ri_,t), where r±_ — xx -\- yy is the two- 
dimensional position vector in the plane perpendicular 
to the z-axis. The interface temperature is then simply 
given by, Tj = T{fint,t), where fint = r± + C{f±,t)z. We 
can assume, without loss of generality, that the interface 
is locally single valued (i.e. no overhang) with respect to 
this set of coordinates; ri(f±,t) is then Gaussianly dis- 
tributed with a variance defined by 



<77(rl,i)77(r'^,t')> = 2 



keTf 6{r^-r'^)5{t-t') 



fiL 



1 + |ViC(^i,t)|^ 



(6) 



where the square-root in the denominator of Eqn. o is 
a simple geometrical factor introduced such that the net 



force on a small area dS of the interface is independent 
of its local orientation |g^, and Vj^ = xd^ +ydy denotes 
the two-dimensional gradient vector in the plane of the 
interface. The application of the fluctuation-dissipation 
theorem for this noise requires that its variance be cho- 
sen such that the sharp-interface model reproduces the 
known fluctuation spectrum of a stationary interface in 
thermal equilibrium, derived analytically in the next sec- 
tion (see also p7|): 



s{k) =< CfcC-fc >= 



7 fc2 



(7) 



where C^j is the Fourier coefficient of the interface dis- 
placement, i.e. 



C(n 



d^k 



Jk-r^ 



Ck 



(8) 



A straightforward but lengthy calculation described in 
Ref. [E^ shows that Eqs. y^, with the noises defined by 
Eqs. Qand ||, yields this spectrum in equilibrium. 



III. PHASE-FIELD MODEL 

The Langevin formalism to incorporate fluctuations 
into continuum models of phase transitions is well- 
established [E6l , and the same procedure can be followed 
for the phase-field model. Like in the sharp-interface 
model Q , we proceed by adding stochastic forces whose 
magnitudes are determined by making contact with equi- 
librium properties. For this purpose, it is convenient to 
express the phase-field model in terms of the dimension- 
less temperature field 



T-T, 



M 



L/c ' 
and the local enthalpy per unit volume defined by. 



H 



Co u 



P(0) 



(9) 



(10) 



where eo is a constant with units of energy per unit vol- 
ume, 4> is the phase- field chosen to vary between —1 
in the liquid and +1 in the solid, and p{(j)) is some 
monotonously increasing function of (f) with the limiting 
values, p(±l) = ±1. The phase-field model expressed in 
terms of these variables takes the form 






fii) 



fl2) 



which is a form similar to Model C of Halperin, Hohen- 
berg, and Ma |^, i.e. with coupled non-conserved (0) 



and conserved (H) order parameters, which is most nat- 
urally suited to add fluctuations. The fact that H is con- 
served, which follows from Eq. 12, simply reflects the fact 



that the total energy in a given volume is conserved in the 
absence of energy fluxes through the surfaces bounding 
this volume. Next, the free-energy is defined by 



Now comparing Eq. n9 with Eq. o allows us to determine 

(20) 



eo 



Tmc 



T = d 



'> 



7/2"! 
P + /,o/(0) + eoy 



(13) 



This result can be obtained, alternatively, by compar- 
ing the phase- field equations (12) and (p^), in a region 
where (p is constant, with the sharp-interface equations 
(|l|) and (jj), which yields, in addition, the expression for 
the diffusion constant 



where Hq and K are constants with units of energy per 
unit volume and per unit length, respectively, and /(</>) is 
a double well potential with minima at = ±1. Specific 
choices for p{cj)) and f{(j)) will be given in the next section 
to carry out the computations. Finally, the noises are 
Gaussianly distributed with variances 



D = 



i_H_ 



(21) 



<e{r,t)9{r',t) > 
< qmir,t)qn{r',t) > 



2T^kBTM5{r- r^)S{t - t') (14) 

2THkBTM 5mn6{r- r')5{t - t') (15) 



Let us now briefly analyze the equilibrium bulk and 
interface fluctuations in this diffuse interface model in 
order to make contact with the sharp-interface model 
of the previous section. As is well-known, the proba- 
bility, P[(/), m; i], of flnding the system in a given configu- 
ration, (/i(r, t) and u = u{f, t), at time t is governed by a 
generalized Fokker-Planck equation |2^ associated with 
the Langevin equations y_l| and |lj. For a general non- 
equilibrium situation, this Fokker-Planck equation has 
no known analytical solution. In equilibrium, however, it 
has a time-independent stationary solution 



P.. 



1 



t>,u\ = -exp 



T 



(16) 



which allows us to calculate analytically the equilibrium 
Gaussian fiuctuations. Here, 



Z = / V(l)Vu exp - 



T 



kuTM 



(17) 



is the equilibrium partition function where !?(/) and "Du 
denote functional integration over the fields (j) and u, re- 
spectively. Let us first calculate the temperature fluc- 
tuations in the bulk phases. Since cj) is constant in the 
solid or liquid, only the term ~ u^ in the integrand of J- 
needs to be kept. Consequently, Eq. Hq implies that the 
fluctuation of u inside a small volume Ay is given by 



< u > = 



du u exp 
dwexp 



AVep u^ 
ksTM T 
AV^eo u^' 
ksTM T 



/ 



(18) 



which yields at once the result, 

ksTm 



< M >: 



eoAl/ 



(19) 



Next, the equilibrium fluctuations of a stationary in- 
terface can be calculated provided that we restrict our 
attention to wavelengths that are large compared to the 
width of the spatially diffuse interface region. Let us con- 
sider the fluctuations about a flat interface in the plane 
z — Q. For a small amplitude deformation, C,{f±^), which 
varies slowly on the scale of the interface thickness, the 
phase-field can be approximated in the form 



(f)«(/.o(z-C(ri)), 



(22) 



where (l)o{z) is the solution of the one-dimensional sta- 
tionary interface problem 



K 



dz^ 



+ hQU{Mz)) = o, 



(23) 



where we have defined /^ = df/d(t). We can then evaluate 
the gradient term in Eq. O using Eq. 22 , which yields 



V0(f) « ^ 
dz 



z^^±C{f±) 



(24) 



Substituting Eq. e3 into Eq. O, we obtain at once that 
the probability distribution of interface fluctuations is 
given by 



p[ar±)] 



1 



■ exp 



7 



ksTM 



dV-|v^C(r^)p 



where 

Z = / PC exp ( - 
and 



ksTM 



d'r-\W^C{r±)\'- 



(25) 



(26) 



7 



y/Kho 



-t-oo 



dz 



dz 



VKh)I (27) 



is the surface energy; the integral / defined in terms of 
the dimensionless variable z = zyJu^jK is a numeri- 
cal constant that depends on the form of /((/>). The re- 
sult of Eq. M stated earlier is now simply obtained by 
changing variables from C(ri) to Cfc in the probability 



distribution above, and by using this distribution to cal- 
culate < CfcC-fe >i which only involves a Gaussian inte- 
gral. This simple exercise shows that the interface fluc- 
tuations in the phase-field model are identical to those 
of the sharp-interface model on scales larger than the 
interface thickness, i.e. k^^ ^ y/K/ho, as one would 
naively expect. Another way to see this correspondence 
is to note that ( p5| ) is identical to the probability distribu- 
tion of small amplitude fluctuations in the sharp-interface 
model. In this model, the probability of an arbitrary fluc- 
tuation of the interface is ~ exp(— 7 J da/fc^TM) where 

da ~ cPrJl -|- |V_lCP is the element of surface area. For 
a small amplitude fluctuation, the square root can be 
expanded to first order which yields at once the distribu- 
tion (p5[). Finally, Eqs. GG and 27 can be used to relate 



the parameters of the phase-field model to the capillary 
length, do = ^TmcJL^ , which yields 



^KU 



(28) 



IV. NUMERICAL IMPLEMENTATION 

A. Choice of functions and scalings 

To carry out numerical simulations, it is convenient to 
choose the functions 






-0V2 + <^V4, 

15((/) - 2(/.V3 + 0V5)/8, 



(29) 
(30) 



where (E9|) is the standard quartic form of double well 
potential and the form (pfl) has the advantage that it 
preserves the minima of at ±1 independently of the lo- 
cal value of u p9|]. The one-dimensional stationary profile 
solution of Eq. |23| is then given by 



00 (-z) = — tanh 



^/2W 



where, 



W 



ho' 



(31) 



(32) 



is the interface thickness. Evaluating the integral in Eq. 
P?! with the above form of (/)o(2:) yields / — 2\/2/3. 

It is useful to express the phase-field equations in a di- 
mensionless form that minimizes the number of computa- 
tional parameters and that renders the interpretation of 
the noise magnitude in the phase-field model more trans- 
parent. For this purpose, it is useful to define, in addition 
to W, the time 



1 



r<A/io 



(33) 



which characterizes the relaxation of < 
minima, and the coupling constant |3 



to one of its local 



A = 



Jho 



I 1 

J Jo' 



expressed in terms of the scaled capillary length 
7 do 



(34) 



(35) 



Here, J — 16/15 is a constant whose value is fixed by 
the choice of p{<j)) ||21|] . We then measure all lengths in 
units of W and time in units of r, and define accordingly 
new dimensionless coordinates, diffusivity, and noise vari- 
ables, via the substitutions 



r/W 

t/r 

Dt/W^ 

T0 

T 



eoW 



t 

D 

9 

q 



(36) 

(37) 
(38) 
(39) 
(40) 



Transforming the phase-field equations O and U2 with 
the help of these substitutions, and using the fact that 
S{f — r') and 6{t — t') on the r.h.s. of equations (nj) 
and (nS) have dimensions of (length)"'^, where d is the 
dimension, and inverse time, respectively, we obtain the 
dimensionless form 



dt 



du 
'dt 

with. 



^ ^DV^u -h - 



•3 -Aw (1-^2)2 ^ g,^-^) (-4-^<) 

1 dp{cj,) 



2 dt 



V-qif,t) 



(42) 



<e{f,t)0{r',t)> ^2F^d{r-r')d{t-t') (43) 

<9™(r,i)9n(r',i)> =2DF^S,-anSif-r^)Sit-t') (44) 



and the definitions, 
p — 

^ ex — 



L^di 



F^^KJR, 



(45) 

(46) 
(47) 



The above definitions allow us to relate the magnitude of 
the noise which enters into the phase-field model, i^^, 
with the magnitude of the noise in experiments, Ff,x- 
Comparing the r.h.s of Eq. Ba, for d = 3, with the 
r.h.s. of Eq. HS, we can readily see that F^^ is simply 
equal to the mean-square fluctuation of u inside a mi- 
croscopic volume do, and is a fixed quantity for a given 
material. (Note that F^x can also be written in the form 
^bTm lld'^, which is the square of the ratio of two micro- 
scopic lengths, y/ksTM M, and do-) The first equality in 



Eq. ^ implies that Fu is the mean-square fluctuation of 
u inside a microscopic volume W^. The second equality 
dictates how to choose F„ in a simulation for a given sys- 
tem (Fex) and a given choice of computational param- 
eter, do. The dependence on the latter quantity has a 
simple physical interpretation. Namely, if one chooses do 
to be small compared to unity, which is the main gain in 
computational efficiency resulting from the reformulated 
asymptotics of Ref. |2^| , then one must scale down the 
magnitude of the noise in the phase-field model to keep 
the fluctuation strength in a physical volume dg constant. 
The main practical conclusion here is that one still has 
the computational freedom to choose the interface thick- 
ness if one rescales appropriately the noise strength. 



B. Discretization 

We first discuss the numerical implementation of the 
model in two dimensions and then briefly mention its 
straightforward extension to three dimensions. The 
phase-field equations mW) and (E2|) are discretized on a 
N X N square lattice of spacing Ax = Az using centered 
finite difference formulae, as described in Ref. [pTl, and 
the equations are time-stepped using a first order Euler 
scheme with a time step At. The only new elements here 
are the noises. To see how to discretize them, let lAx 
and jAz denote the position on the lattice along x and 
z, respectively. For the non-conserved noise, we generate 
one random number per lattice site, 0^, chosen from a 
Gaussian distribution with a variance 

2F^ 



< 



'V^i'j' 



> = 



AtAa 



^ii' ^jj' 1 



(48) 



where the factors 1/ At and 1/Aa;^ on the r.h.s. of Eq. 
Bq are related to the inverse time and the inverse area 
(inverse volume in three dimensions) scalings oi 5{t — t') 
and 5{r — r'), respectively, in the correlation of the noises. 
9ij is then added to the deterministic part of the r.h.s. 
of Eq. El] discretized at site («,j). 

To discretize the conserved noise, we define by qx,ij 
the current on the bond that links site {i,j) with site 
{i + 1, j), and by qz^ij the current on the bond that links 
site (i,j) and (i,j -t- 1). We then generate at each time 
step two independent random numbers per site, qx,i] and 



qz 



chosen from a Gaussian distribution with a variance 



<* 



771 ^IJ Hn,i J 



> 



2DFu 



^mn^ii' ^i 



AtAx"^ 

The divergence of the current on the r.h.s. of Eq. 
site (i, j) is then discretized in the form 



(49) 



laat 



[V ■ qj = [qx,ij - q. 



'x,i-lj 



dz,tj 



qz,^J-l]/Ax (50) 



In three-dimensions, the only changes involve replacing 
Ax^ by Ax^ in the denominator of the r.h.s. of equations 
(El) and (|49|), and to generahze (pO[) to 



I V • (f 1 — [qx.ijk — qx,i-ljk + qySjk — qy,ij-lk 

\ / ijk 

+qz,ijk - qz.i]k-i] /Ax (51) 



C. Planar interface fluctuation spectrum 

As a non-trivial test of the numerical implementation 
of the phase-field model, we first calculate the fluctua- 
tion spectrum of a one-dimensional stationary interface 
in thermal equilibrium and compare this spectrum to the 
analytical prediction (|^). With length measured in units 
of W, (0) becomes 



S{k) = 



dok'^ 



(52) 



where Jo — I /{J^) = 5\/2/8A for the present choice of 
phase-field model. 

To calculate S'(fc), phase-field simulations were carried 
out with periodic boundary conditions in x on a lattice of 
size 512 X 50 with Ax — 0.8. We used initial conditions 
that correspond to a flat interface inside a system uni- 
formly at the melting temperature, which corresponds to 
choosing — (l)o{z) = — tanh(z/\/2) and m = 0. The 
interface profile, C,{x)^ is defined by (j){x,C,{x)) = 0, and 
is calculated by finding the (f) — Q contour of the phase- 
field by interpolation at each time step. The complex 
amplitude, Cfc, is then calculated by a one-dimensional 
fast Fourier transform where Qk and Q{x) are related by 



C(^) = / § e^"' 



Ck, 



(53) 



Finally, S{k) —< |Cfep > is calculated by taking a time 
average of |CfcP- Long simulations with typically 10^ to 
10^ time steps were necessary to obtain good statistics. 
These calculations were carried out by using (E2|) with 
both p{(j}) defined by (|^) and p(0) = cj). With the latter 
choice, the phase-field equations arc no longer variational 
(i.e. derivable from a Lyapounov functional), but the 
sharp-interface limit remains identical and the interface 
can be resolved with a larger Ax, as shown previously 
pi| . The spectra for the two choices of p((?!)) were found 
to be virtually indistinguishable such that only the re- 
sults for p{4>) = 4> are reported here. In the dendritic 
growth simulations presented below, we will restrict our 
attention solely to the case where p(0) = is used as the 
source of latent heat in the heat equation. 

Spectra obtained for a typical set of computational pa- 
rameters are compared in Fig. || with the analytical pre- 
diction (p2|), represented by a thick solid line, both for 
the case where the non-conserved and conserved noises 
are added to the phase-field equations (thin solid line 
with Fu y^ and F^ ^ 0) and for the case where the non- 
conserved noise is switched off in the (j) equation (dashed 
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FIG. 1. Simulated spectra of a one-dimensional interface in 
thermal equilibrium with both non-conserved and conserved 
noises (thin solid line) and only the conserved noise (dashed 
line) , compared to the theoretical prediction of Eq. B2l (thick 
solid line). Length is measured in units of W. Parameters 
used in simulations are A = 1, D = 1, and Fu — 0.005. 



line with F^ ^ but F^ = 0). With both noises present, 
the calculated spectrum agrees well with the theoretical 
prediction up to a cutoff in k of order unity (correspond- 
ing to a wavelength comparable to the interface thickness 
in physical units). With only the conserved noise present 
{F^ = 0), the simulated spectrum follows initially well 
the predicted 1/fc^ law with increasing fc, but then drops 
off rapidly to a very small amplitude at large k. This drop 
off is consistent with the analytical prediction of Ref. ||27| ] 
and is due to the extra dissipation at the interface that 
damps out short scale fluctuations. 



D. Incorporation of anisotropy 

In order to investigate dendritic sidebranching in the 
next section, we incorporate anisotropy as other authors 
pl| , p2[ by letting the coefficient of the gradient energy 
term in the free-energy depend on the normal to the solid- 
liquid interface, fi — V0/|V0|. Following this change, 
(E2) remains unchanged and ([4l|) becomes 



2\2 



/fe(n)9t0 = 0-0^ -A 7/(1-0^) 



^V-(/,(n)2V0) 
e{r,t) 



m—x,z 



dm \^<jy?fs{n: 



dfsjn) 
d{dm(f>) 



(54) 



where we have defined the anisotropy function for a crys- 
tal with an underlying cubic symmetry 



fs{n) = 1 - 364 + ie^idU + dU)/\^cl)\' 



(55) 



As in our previous study of dendritic growth without 
noise pll] , we use the result of a reformulated asymptotic 



TABLE L List of the phase-field computational parameters 
used in dendritic growth simulations. These parameters yield 
an effective 3% anisotropy in surface energy and a diverging 
interface kinetic coefficient fi as defined here in Eq. H. 



Aa; 

At 
D 
A 
do 

£4 



0.8 
0.06 

2 
3.268 
0.27 
0.03266 
0.03 
0.046 



analysis of the phase-field model together with a method 
to compute lattice corrections to the surface energy and 
kinetic anisotropies. Moreover, we focus on a choice of 
computational parameters that makes 1/ ^ vanish in the 
interface condition (||). The effective anisotropy of the 
phase- field model, which includes lattice corrections, is 
given at order Ax^ by 



Ce = £4 - A2;V240 



(56) 



Here, we use Ax = 0.8, and input the value £4 = 0.03266 
into Eq. pq to obtain an effective 3% anisotropy when 
comparing our results to the sharp-interface solvability 
theory. To make the interface kinetic contribution vanish 
at the same order we choose 



fkin) = 



l-35 + 4S{dt<l> + dt(l))/\W<j>\' 
l + S 



(57) 



where the value of S is computed, together with an or- 
der Ax'^ correction to A, in order to make l//i vanish 
in (y|), as described in |gl|]. The resulting computational 
parameters are summarized in Table I. 

Lastly, in terms of our dimensionless units, where 
length, time, and velocity, are scaled in units of W , r, 
and W/t, respectively, and without interface kinetics, 
the thin-interface limit of the phase-field model is the 
standard free-boundary problem: 

dtu = DV^u-V -q (58) 

Vn = -Dh- (Vu|, - Vu|,) + h-{q]i^q]s) (59) 

uj = ~do(l ^ ISce COS 4a) k (60) 

where g is the same noise as in the phase- field model and 
a = cos~^ (z-fi) is the angle of the normal measured from 
the z-axis. 



V. DENDRITIC SIDEBRANCHING 

In this section, we simulate the phase-field model de- 
fined by Eqs. E3 and M to investigate sidebranching 



characteristics for different noise levels and a fixed di- 
mensionless undercooling A = {Tm — Toa)/{L/c) = 0.55, 
where Too is the initial temperature of the melt. We then 
compare these results quantitatively with the predictions 
of the linear WKB theory that corresponds to the sharp- 
interface model defined by Eqs. 



A. Numerical results 

Test simulations were first carried out with both noises, 
6 and q, and with only the conserved noise q. We 
found that time-averaged sidebranching characteristics 
were identical for the two cases within our numerical 
resolution. This finding shows that fluctuations which 
become amplified to produce sidebranches are on length- 
scales much larger than the interface thickness, and thus 
driven solely by the bulk noise in agreement with ex- 
pectation (see section I and ^^). Consequently, all the 
results presented in this section were obtained with sim- 
ulations where noise is added only to the heat transport 
equation (p2). This represents a non-negligible compu- 
tational saving for long simulation runs (i.e. 2 instead of 
3 random numbers per site at each time step). 

The development of a dendrite and its sidebranches 
from a small initial seed is illustrated in Figs. ||a and 
pp for two different noise levels. These particular sim- 
ulations were carried out on a large 1200 x 1200 lattice 
with, as initial condition, u — Q and 0=1 inside a 
small circle in the lower left hand corner of the quad- 
rant and u — —A and ip — —1 outside this circle. Note 
that in Fig. pb the noise is sufficiently large to disturb 
the steady-state growth of the tip, which can be deduced 
from the fact that the vertical branch has outgrown the 
horizontal branch in this case. Since the tips do not in- 
teract via the diffusion field at this undercooling, i.e. the 
separation between the tips is much larger than the dif- 
fusion length, this difference can only be due to noise. 
This effect is negligible for the smaller noise level (Fig. 
ya) where the two tips grew at nearly the same rate. 

To investigate sidebranching, we restrict our atten- 
tion to small noise levels {F^ — 2.5 x 10"^ and F„ = 
2.5 X 10~^) with a well-defined steady-state tip struc- 
ture. The symmetric growth of one tip about the z-axis 
(i.e. half the dendrite with reflection symmetry) is simu- 
lated on lattices of size 300 x 400 and 300 x 600 for, re- 
spectively, the larger and smaller noise amplitude (where 
sidebranches form further behind the tip). As in Ref. 
PH , we periodically translate the entire structure in the 
opposite direction of growth to allow long simulation runs 
to be carried out in the smallest lattice size possible. Of 
course, we make sure that the sidebranching activity is 
not affected by this procedure by choosing a reasonable 
buffer larger than the diffusion length, and by carrying 
out test runs with larger lattice sizes. The constraint 
of symmetric growth prevents us from investigating the 
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FIG. 2. Morphological development of a solid seed for 
A = 0.55 and two different conserved noise amplitudes: (a) 
Fu — 10~*, and (b) Fu = 10""^. Other computational param- 
eters are listed in Table 1. The interfaces are plotted every 
16000 iterations. 



correlation of the sidebranching activity on opposite sides 
of the growth axis, which has been examined experimen- 
tally. However, it permits a more efficient investigation of 
the sidebranching amplitude and wavelength which can 
be compared to analytical predictions. 

To calculate these two quantities, we proceed in two 
steps. Firstly, we carry out a simulation without noise to 
obtain a 'reference' steady-state (needle crystal) shape, 
without sidebranching. It is useful to measure this shape 
by the horizontal distance, xo{z), of the interface mea- 
sured from the vertical growth axis as a function of the , 
distance, z, behind the tip. This distance is calculated 
by numerical interpolation of the (f> = contour. The 
steady-state operating state of the dendrite is defined in 
terms of the dimensionless tip velocity and radius 



V = Vdo/D 
P = p/do 



(61) 
(62) 



For the present choice of undercooling and anisotropy, we 
find that V = Vdo/D ^ 0.011 and p = p/do « 21.8, in 
good agreement with the predictions of solvability the- 
ory 1^,^ ^. Secondly, we add the conserved noise 
to the heat equation and calculate the time-dependent 
shape, x{z,t), with sidebranching present. Snapshots of 
noisy shapes superimposed on the noiseless shape are il- 
lustrated in Fig. 0. Examples of time traces of x{z,t) 
for two different distances behind the tip are shown in 
Fig. H. In addition, an example of the noise averaged 
power spectrum of a long time trace is shown in Fig. |g. 
This spectrum was calculated by subdividing the com- 
plete time interval into several equal subintervals, then 
calculating the power spectrum for each subinterval, and 
finally taking the average of these power spectra. 

In terms of the above quantities, the root-mean-square 
amplitude of sidebranches is simply given by: 



A{z) = v/< ixiz,t)^xoiz)y > 



(63) 



where the average is over time. This quantity is plotted 
vs z in Fig. for two different noise levels. To obtain 
good statistics, we typically simulated a total time of 
2000V/ p which took 200-350 CPU hours on a high end 
workstation. The mean spacing between sidebranches 
(sidebranching wavelength), < A(z) >, can be calculated 
in two ways. One way, which corresponds more directly 
to the way in which this quantity is calculated in the 
WKB theory discussed below, is to define 



< A(z) > = 



2'kV 

LUc{z) 



(64) 



where lUc is the peak frequency of the power spectrum of 
x{z, t), averaged over sufficiently long time. An alternate, 
and simpler way, which avoids to calculate the power 
spectrum, is to count the number, N{z), of extrema of 
x{z,t) in a long time interval ti < t < 12- Simple node 
counting then leads to the relation 
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FIG. 3. Snapshots of the time-dependent dendrite shapes 
(solid lines) in long simulation runs that focus on the growth 
of one tip for A = 0.55 and the parameters of Table I. The 
noiseless shape (dashed line) is superimposed for comparison. 
The noise levels are Fu = 2.5x10"^ in (a) and F„ = 2.5x10"'' 
in (b). Note that sidebranches form further behind the tip for 
the smaller noise level. 
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FIG. 4. Horizontal position of the interface measured from 
the vertical growth axis as a function of dimensionless time 
at 20 and at 40 tip radii behind the tip. The parameters are 
the same as in Fig. Hb. 
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FIG. 5. Noise averaged powerspectrum of x{z, t) at 
|2|/p = 40 for F^ = 2.5 x 10"^ 
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FIG. 6. Root-mean-square amplitude of sidebranches as 
a function of distance behind the tip for two different noise 
levels in the simulations (thick solid and dash lines). Super- 
imposed are thin solid lines corresponding to the analytical 
predictions of Eq. 66l 



< \{z) > = 



N{z) 



(65) 



We have checked that these two methods yield nearly 
identical values for the spectrum of Fig. Q Therefore, 
we have used Eq. |6^ to calculate < X{z) > vs z and the 
result for the lowest noise level is shown in Fig. Q 



B. Comparison with linear WKB theory 

Langer ]10| , and Brener and Temkin [ pdj , have an- 
alyzed noise-induced sidebranching in three dimensions 
for specific needle crystal shapes (i.e. x ~ z^/^ and 
X ~ z^^^). It is straightforward to extend their analyses, 
based on a WKB approach, to an arbitrary needle crys- 
tal shape, xo{z), in d-dimension |g^. We shall only state 
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FIG. 7. Mean spacing of sidebranches as a function of dis- 
tance behind the tip in the simulation for Fu — 2.5 x 10~® 
(solid line), analytically predicted by Eq. pi (dashed line) 
and predicted on the basis of stretching (dotted line). 



here the final results necessary to interpret our simula- 
tions. The expressions for the sidebranching amplitude 
and wavelength are given, respectively, by p7[|: 



A{z) = S* exp I - 



3a*z 



1/2N 



< A(z) > = ^ 



12a* 



-,1/2 



XO 



(66) 



(67) 



where we have defined the scaled quantities xq = xo/p, 
z = —zjp, A = A/ p, A = A/p, u) = ujp/V , the dimen- 
sionless noise amplitude, S*, given by 



S' 



2F,, D 



2F„ 



and 



P^+'^V d^pT^+dy 



2Ddo 
pW 



(d = 2,3) (68) 



(69) 



It is easy to check that Eqs. pO and p^ reduce to the 
earlier results of Refs. po| , pl| if specific shapes (parabola 
and 3/5 law) are substituted into them. Note that if we 
convert back to dimensional units by letting p —^ p/W , 
V -^ Vt/W, and D -^ Dt/W'^, in Eq. ||, we obtain 
the expression 



S' = 



2kBTljcD 
L^p^+dy 



(70) 



which is dimensionless if one interprets L and cTj^i to 
have dimension of energy/ (length)''. Of course, this in- 
terpretation is only physically meaningful in three dimen- 
sions where Eq. pfl becomes identical to the definition 
of S in Ref. pOl. Therefore, in the present study we 
evaluate S directly from Eq. p8^ to compare simulations 
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and theory. For do = 0.27 (table I) and the aforemen- 
tioned selected values, V « 0.011 and p « 21.8, we ob- 
tain S « 0.24 F„ and a* « 0.383. (Note that a* is larger 
here than in experiment due to both, the large value of 
anisotropy which produces a pointy tip, and the fact that 
a* is larger in 2-d than 3-d for the same anisotropy.) 
The analytical predictions for the sidebranching ampli- 
tude and wavelength are then simply obtained by in- 
puting these values into Eqs. ^ and ^ together with 
the steady-state interface shape, xo{z), measured in the 
noiseless phase-field simulation (dashed lines in Fig. ||). 

Figs. ^ and ^ show that the amplitude and wavelength 
measured in the phase-field simulations with noise are in 
good overall quantitative agreement with the analytical 
predictions even though a* is not much smaller than one. 
The amplitude in the simulations is relatively well pre- 
dicted by Eq. |6^, up to a certain distance behind the tip 
after which the two curves depart from each other. This 
departure may be due to nonlinear effects which become 
important when A becomes of order unity. In addition, 
it should be emphasized that the prefactor of Eq. ^ 
is only known up to some multiplicative factor of order 
unity. Consequently, what is more relevant here is that 
the amplitudes in simulation and theory are of compara- 
ble magnitude, rather than the fact that the numerical 
and theoretical curves in Fig. g seem to almost perfectly 
overlap up to z ~ 20, which may be coincidental. 

The wavelength in the simulations is only about 30% 
larger than predicted by Eq. |6^ in the region (20 to 
40 tip radii behind the tip) where sidebranches becomes 
visible. However this wavelength increases initially faster 
with distance behind the tip than predicted. One possible 
explanation for this faster rate of increase is that pertur- 
bations generally get stretched as they travel along the 
sides of curved fronts @-^ . To test this possibility, let us 
calculate the purely deterministic change of wavelength 
of a perturbation initially at the tip due to stretching. 
The rate of stretching is given by 0-0] : 



IdX 

xH 



ds 



(71) 



where Vt — V sin a is the tangential velocity of advec- 
tion of the perturbation and s measures the arclength 
along the interface. Eq. ^ is strictly valid in the WKB 
limit where A is small compared to the local radius of 
curvature (1/k) of the interface. We can solve Eq. |7l| 
by using the change of variable dt = dz/V. Eq. ^ be- 
comes then, d(lnA) — sin2adQ;/2, which can be easily 
integrated. Furthermore, using the geometrical relation, 
cos a = 1/[1 + (dxo/dz)^]"^/^, we obtain 



Xiz) = Aoo exp 



1 [dxp/dzf 

2 l + {dxo/dzy 



(72) 



where we have defined Aoo to be the saturation value of 
A far behind the tip (i.e z — > cxd). In order to test if the 



disagreement between simulation and theory in Fig. 
is due to stretching near the tip, we have plotted in the 
same figure the prediction of Eq. [T^ with Xoc/p = 13.5 
fitted from the simulation data. We can see that this pre- 
diction indeed improves the agreement with simulation in 
this region. This stretching effect is not contained in Eq. 
pT) which is only valid in the far tip region {z ^ 1). The 
slower increase in < A > predicted by Eq. ^is due to the 
distinct effect that lower frequency perturbations survive 
at a larger distance from the tip &M. It should be in 
principle possible to carry out a more elaborate WKB 
analysis with noise that includes both the latter effect 
and stretching. Such an analysis will yield the same pre- 
diction as Eq. |67| in the far tip region, and (presumably) 
an improved wavelength prediction closer to the tip. 

Finally, we note that the sidebranching wavelength is 
about an order of magnitude larger than the tip radius in 
our simulations, whereas it is only a factor of two or three 
in the classic experiments of Huang and Glicksman in 
sunccininitrile 0. This difference is due to the fact that 
a* is much larger here than in these experiments because 
of the larger value of anisotropy used in simulations. 



VI. CONCLUSIONS 

We have presented a phase-field model of the solidi- 
fication of a pure melt that incorporates thermal noise 
quantitatively. From a computational standpoint, there 
are two main conclusions regarding the incorporation of 
this noise. Firstly, one can retain the freedom to choose 
the interface thickness at will as long as the noise mag- 
nitude (Fu) that enters in the phase- field model is scaled 
appropriately (Eq. Em. Therefore, it remains possible 
to carry out dramatically more efficient computations 
without interface kinetics by choosing do substantially 
smaller than unity, as in our earlier studies without noise 
[|2l| . Secondly, for typical growth conditions at low un- 
dercooling (and, more generally, below a critical velocity 
that depends on the attachment kinetic coefficient /i), the 
conserved noise in the heat current is the most relevant 
one. This noise drives long-wavelength interface fluctu- 
ations that become amplified to a macroscopic scale by 
the morphological instability on the sides of dendrites. In 
contrast, the non-conserved noise in the evolution equa- 
tion for (j) drives short-wavelength fluctuations that are 
damped and do not affect sidebranching, as predicted 
by a sharp-interface analysis [ p7| and confirmed by our 
simulations. Consequently, this noise can be left out in 
computations below this critical velocity. 

We have applied this model to carry out a detailed 
quantitative study of the initial stage of sidebranch for- 
mation during dendritic growth. The main conclusion is 
that the sidebranching characteristics (root-mean-square 
amplitude and the mean sidebranch spacing) are reason- 
ably well predicted by the existing linear WKB theory 
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of noise amplification, even though the value of a* in 
our simulations is not small. A more stringent test of 
this theory would therefore require to extend the present 
study to a range of smaller anisotropy, and hence smaller 
(J* , where a comparison between WKB theory and sim- 
ulations is more justified. Finally, since this study has 
been restricted to two dimensions, we cannot yet answer 
the important question of whether thermal noise alone 
is responsible for the sidebranching activity observed in 
experiment. Simulations in three dimensions should pro- 
vide a clear cut answer to this question. 
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